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Using computer simulations we investigate the microscopic 
structure of the singular director field within a nematic 
droplet. As a theoretical model for nematic liquid crystals we 
take hard spherocylinders. To induce an overall topological 
charge, the particles are either confined to a two-dimensional 
circular cavity with homeotropic boundary or to the surface of 
a three-dimensional sphere. Both systems exhibit half-integer 
topological point defects. The isotropic defect core has a ra- 
dius of the order of one particle length and is surrounded 
by free-standing density oscillations. The effective interac- 
tion between two defects is investigated. All results should 
be experimentally observable in thin sheets of colloidal liquid 
crystals. 

PACS: 61.30.Jf, 83.70.Jr, 77.84.Nh 



I. INTRODUCTION 



Liquid crystals (LC) show behavior intermediate between 
liquid and solid. The coupling between orientational and 
positional degrees of freedom leads to a large variety of 
mesophases. The microscopic origin lies in anisotropic 
particle shapes and anisotropic interactions between the 
particles that constitute the material. The simplest, most 
liquid-like LC phase is the nematic phase where the par- 
ticles are aligned along a preferred direction while their 
spatial positions are, like in an ordinary liquid, homo- 
geneously distributed in space. The preferred direction, 
called the nematic director, can be macroscopically ob- 
served by illuminating a nematic sample between crossed 
polarizers. 

There are many different systems that possess a nematic 
phase. Basically, one can distinguish between molecular 
LCs where the constituents are molecules and colloidal 
LCs containing mesoscopic particles, e.g., suspensions of 
tobacco mosaic viruses [|l| . Furthermore there is the pos- 
sibility of self-assembling rodlike micelles Q , that can be 
studied with small-angle neutron scattering [||. 

There are various theoretical approaches to deal with ne- 
matic liquid crystals. On a coarse-grained level one may 
use Ginzburg-Landau theories, including phenomenolog- 
ical elastic constants. The central idea is to minimize 



an appropriate Frank elastic energy with respect to the 
nematic director field [Q. Second, there are spin mod- 
els, like the Lebwohl-Lashcr model, see, e.g., Refs. [pj-R. 
There the basic degrees of freedom are rotators sitting 
on the sites of a lattice and interacting with their neigh- 
bors. The task is to sample appropriately the configu- 
ration space. The third class of models consists of par- 
ticles with orientational and positional degrees of free- 
dom. Usually, the interaction between particles is mod- 
elled by an anisotropic pair potential. Examples are Gay- 
Berne particles, e.g. [ffljl and hard bodies, e.g., hard 
spherocylinders (HSG) |lO||. Beginning with the clas- 
sical isotropic-nematic phase transition for the limit of 
thin, long needles due to Onsager ||ll| , our knowledge has 
grown enormously for the system of hard spherocylinders. 
The bulk properties have recently been understood up to 
close packing. The phase diagram has been calculated 
by computer simulations |12|, density- functional theory 
p3[ and cell theory flj]. There are various stable crys- 
tal phases, like an elongated face-centered cubic lattice 
with ABC stacking sequence, a plastic crystal, smectic-A 
phase, nematic and isotropic fluid. Besides bulk proper- 
ties, one has investigated various situations of external 
confinement, like nematics confined to a cylindrical cav- 
ity ||l5| or between parallel plates |16 1^. Also effects 
induced by a single wall have been studied, like depletion- 
driven adsorption ||l^], anchoring [^, wetting [^, and 
the influence of curvature pl|. Furthermore, solid bodies 
immersed in nematic phases experience non-trivial forces 
22I-M , and point defects experience an interaction pSl . 



Topological defects within ordered media are deviations 
from ideal order, loosely speaking, that can be felt at an 
arbitrary large separation distance from the defect po- 
sition. Complicated examples are screw dislocations in 
crystalline lattices and inclusions in smectic films [p6|. 
To deal with topological defects the mathematical too s 
of homotopy theory may be employed [E^ to classify all 
possible structures. The basic ingredients are the topol- 
ogy of both the embedding physical space and the order 
parameter space. For the case of nematics, there are two 
kinds of stable topological defects in 3d, namely point de- 
fects and line defects, whereas in 2d there are only point 
defects. These defects arise when the system is quenched 
from the isotropic to the nematic state pS] . Also the dy- 
namics have been investigated |29| ] experimentally. On 
the theoretical side, there is the important work within 



the framework of Landau theory by Schopohl and Sluckin 
on the defect core structure of half-integer wedge dischna- 
tions Q ^^'^ ^^ ^^^ hedgehog structure |^ in nematic 
and magnetic systems. The latter predictions have been 
confirmed with computer simulations of lattice spin mod- 
els 13^ . The topological theory of defects has been used 
to prove that a uniaxial nematic either melts or exhibits 
a complex biaxial structure |33| ]. Sonnet, Kilian and Hess 
p4| have considered droplet and capillary geometries us- 
ing an alignment tensor description. 

The investigation of equilibrium topological defects in ne- 
matics has received a boost through a striking possibil- 
ity to stabilize defects by imprisoning the nematic phase 
within a spherical droplet. The droplet boundary induces 
a non-trivial effect on the global structure within the 
droplet. Moreover, it can be experimentally controlled in 
a variety of ways to yield different well-defined boundary 
conditions, namely homeotropic or tangential ones. One 
famous experimental system are polymer-dispersed LCs. 
Concerning nematic droplets, there are various studies 
using the Lebwohl-Lasher model [pbR- There are inves- 
tigations of the droplet shape [^,0 , the influence of an 
external fiel d p7|| , and chiral nematic droplets [pq] , struc- 
ture factor |3Sf| 7 and ray propagation ||40(| . Also simula- 
tions of Gay-Berne droplets have been performed pi. 
Other systems that exhibit topological defects are ne- 
matic emulsions p2|-p4|, and defect gels in cholesteric 
LCs [ESJ . The formation of disclination lines near a free 
nematic interface was reported Q]. 

In this work we are concerned with the microscopic struc- 
ture of topological defects in nematics. We use a model 
for rod-like particles with a pair-wise hard core interac- 
tion, namely hard spherocylinders. It accounts for both, 
the oricntational degrees of freedom as well as the posi- 
tional degrees of freedom of the particles constituting the 
nematic. Especially, it allows for mobility of the defect 
positions. This system is investigated with Monte Carlo 
computer simulations. There exist successful simulations 
of topological line defects using hard particles, namely 
integer [^ and half- integer line defects pq ]. 

Here, we undertake a detailed study of the microscopic 
structure of the defect cores focusing on the behavior 
of the local nematic order and on the density field, an 
important quantity that has not been studied in the lit- 
erature yet. As a theoretical prediction, we find that 
the arising half-integer point defects are surrounded by 
an oscillating density inhomogeneity. This can be ver- 
ified in experiments. We also investigate the statistical 
properties of two defects interacting with each other ex- 
tracting the distribution functions of the positions of the 
defect cores and their orientations. These are not acces- 
sible in mean-field calculations. We emphasize that both 
properties, the free-standing density wave which is due 
to microscopic correlations and the defect position dis- 
tribution which is due to fluctuations cannot be accessed 



by a coarse-grained mean-field type calculation. 

The paper is organized as follows: In section |l| our the- 
oretical model is defined, namely hard spherocylinders 
within a planar spherical cavity and on the surface of a 
sphere. For comparison, we also propose a simplified toy 
model of aligned rods. Section HI is devoted to the ana- 
lytical tools employed, such as order parameter and den- 
sity profiles. Section |^ gives details about the computer 
simulation techniques used. The results of our investiga- 
tion are given in section M and we finish with concluding 
remarks and a discussion of the experimental relevance 
of the present work in section VI , 



II. THE MODEL 



A. Hard Spherocylinders 



We consider N identical particles with center-of-mass 
position coordinates r^ = {rxi,ryi) and orientations rij, 
where the index i = 1, . . . ,N labels the particles. Each 
particle has a rod-like shape: It is composed of a cylin- 
der of diameter a and length L ~ a and two hemispheres 
with the same diameter capping the cylinder on its flat 
sides. In three dimensions (3d) this geometric shape is 
called a spherocylinder, see Fig.Q^ . The 2d analog is 
sometimes called discorectangle as it is made of a rect- 
angle and two half discs. We assume a hard core in- 
teraction between any two spherocylinders that forbids 
particle overlap. Formally, we may write 



U{r,,n,;rj,nj) 



cxD if particles i and j overlap , _. ^ 
else ^ ' 



The geometric overlap criterion involves a sequence of el- 
ementary algebraic tests. They are composed of scalar 
and vector products between the distance vector of both 
particles and both orientation vectors. The explicit form 
can be found e.g. in Ref. Q. The bulk system is gov- 
erned by two dimensionless parameters, namely the pack- 
ing fraction 77, which is the ratio of the space filled by the 
particle "material" and the system volume V. In two di- 
mensions it is given by 77 = {N/V){a{L — a) + 7rcr^/4). 
The second parameter is the anisotropy p = L/cr which 
sets the length-to- width ratio. The bulk phase diagram in 
3d was recently mapped out by computer simulation 112] 
and density- functional theory [ |l3[ . The nematic phase 
is found to be stable for anisotropics p > 5. In 2d the 
phase diagram is not known completely but there is an 
isotropic to nematic phase transition for infinitely thin 
needles JS^]- The nematic phase is also present in a sys- 
tem of hard ellipses pi 52 1 verified by computer simula- 



tions. In 2d the ncmatic-isotropic transition was inves- 
tigated using density-functional theory J53| and scaled- 
particle theory pM. There is work about equations of 



state [|5|, and direct correlation functions IQ within a 
geometrical framework. 



B. Planar model 



To align the particles near the system boundary 
homeotropically we apply a suitably chosen external po- 
tential. The particles are confined within a spherical 
cavity representing the droplet shape. The interaction 
of each HSC with the droplet boundary is such that the 
center of mass of each particle is not allowed to leave the 
droplet, see Fig.El The corresponding external potential 
is given by 



C^oxt(rj) 



if \vi\<R-L/2 
oo else 



(2) 



where R is the radius of the droplet and we chose the 
origin of the coordinate system as the droplet center. 
The system volume is, V = ttR^ . This boundary con- 
dition is found to induce a nematic order perpendicular 
to the droplet boundary as the particles try to stick one 
of their ends to the outside Q . Hence the topological 
charge is one. In the limit, p = 1, we recover the confined 
hard sphere system recently investigated in 2d |5q] and 



D. Aligned Rods 



To investigate pure positional effects we study a further 
simplified model where the orientation of each rod is 
uniquely determined by its position. Therefore we con- 
sider an arbitrary unit vector field n(r) describing a given 
nematic order pattern. In reality, the particles fluctuate 
around this mean orientation. Here, however, we neglect 
these fluctuation by imposing n^ = n(ri). In particular, 
we chose the director field to possess a singular defect 
with topological charge t, see Fig.H. The precise defini- 
tion of this director field n*^*^ (r) is postponed to the next 
section (and given therein in Eq||.) The case of parallel 
aligned rods, n = const, has been used to study phase 
transitions to higher ordered hquid crystals p^. 



III. ANALYTICAL TOOLS 



A. Order parameters 



In order to analyze the fluctuating particle positions and 
orientations, we probe against a director field possessing 
a topological defect with charge t. It is given by 



n^'Mq,r)-D^^Ur)q, 



(t), 



(5) 



where the rotation matrix is 



C. Spherical model 



A second possibility to induce an overall topological 
charge is to confine the particles to a non-planar, curved 
space, which we chose to be the surface of a sphere 
in three-dimensional space. The particles are forced to 
lie tangentially on the sphere with radius R, see FigJ|. 
Mathematically, this is expressed as 



■ n, 



R, 
= 0. 



(3) 

(4) 



The director field on the surface of a sphere has to have 
defects. This is known as the "impossibility of comb- 
ing a hedgehog". The total topological charge [^ is 
two. The topological charge is a winding number that 
counts the number of times the nematic director turns 
along a closed path around the defect. It may have pos- 
itive and negative, integer or half-integer values, namely 
0,±1/2,±1,.... 



(t)f.\- 



B}''ia) 



cos (tcf)) 
sin (tcj)) 



— sin (i0) \ 
cos {t4>) J 



(6) 



with (p = aT:cta,n{ay / ttx) , and a = {ux^ay) being a 2d 
vector. The vector q is the orientation of particles if one 
approaches the defect along the a;-direction. 

As an order parameter, we probe the actual particle ori- 
entations rij against the ideal ones 



5(*)(c,q;r)=2(k •nW(q,r,;-c) 



-1, 



(7) 



where the radial average is defined as {■ ■ ■)j. = 

{EliS{\r[\-r)...)/(j:liSi\v[\-r)), with r[ = 

rj — c and (. . .) is an ensemble average. Normalization 
in EqJ^ is such that usually < S**^*-* < 1, where unity 
corresponds to ideal alignment, and zero means complete 
dissimilarity with the defect of charge t at position c and 
vector q, Eq.|[ (In general, — 1 < S'^*^ < 1 is possible, 
where negative values indicate an anti-correlation.) 



If c and q are not dictated by general symmetry con- 
siderations (e.g. c = because of the spherical droplet 



shape), we need to determine both quantities. To that 
end we measure the similarity of an actual particle con- 
figuration compared to a defect, Eq.^. We probe this 
inside a spherical region around c with radius R* using 



I^'Hc,<i) ^ j^ J^ drrS^'\c,<i;r). 



(8) 



Second, we probe for star-Hke order, hence t = 1, c = 0. 
As we do not expect spiral arms of the star pattern to 
occur, we can set q = e^;, where e^ is the unit- vector in 
x-direction. We can rewrite Eq.0 as 



5«(r) = 2((n,.f,)'\ -1, 



(14) 



where R* is a suitably chosen cutoff length. We maxi- 
mize /^*-'(c, q) with respect to c and q. The value at the 
maximium is 



A(*) =max{/(*)(c,q)}, 

c,q 



(9) 



and the argument at the maximum is q' 



(t) 



Before summarizing the quantities we compute during 
the simulation, let us note that q^*^ and A*^*-* are eigen- 
vector and the corresponding (largest) eigenvalue of a 
suitable tensor. To see this, we attribute each particle 
the general tensor 

Q,(*)=2fD(*)(r, -c)n,®D(*)(r, -c)n,) -1, (10) 



where (8) denotes the dyadic product, 1^ is the identity 
matrix. Summing over particles gives 



Q 



(*) 



Eq. 



(t) 



(11) 



Note that for t — the usual bulk nematic order param- 
eter is recoveredR The order parameter profile, Eq.R, is 
then obtained as 



^(*)(c,q,r) = (q.2(*).q 



(12) 



where f^ = ri/\n\. 

Third, we investigate t = 1/2 defects. To that end, we 
need to search for c and q, as these are not dictated by 
the symmetry of the droplet. Hen ce we numerically solve 
Eq.| with R* = 2L (see Sec. |[VB| .) We obtain 



^(1/2) (r) = 2 ( (n, • n(i/2) (q^/s) , r, - c^ 



1/2; 



- 1. 



(15) 



The distribution of the positions of the particles is ana- 
lyzed conveniently using the density profile p{r) around 
c, which we define as 



p(r) = /(27rr)-ilf^5(|r,-c|-r) 



(16) 



We consider two cases: The density profile around the 
center of the droplet, i.e. c = 0, and around the position 
of a half-integer defect, c = Ci, C2. 

It is convenient to introduce a further direction of a i = 
1/2 defect by 



d = D(i)fq(i))q(^ 



(17) 



and then the relation A'^'-'q'-*-' — Q^'-'q'*' holds, if the sum 
over i in Eq, 11 is restricted to particles located inside a 



spherical region of radius R* around c. 

Let us next give three combinations of t, c, q that apply 
to the current model. First, we investigate the (bulk) 
nematic order, t — 0. We resolve this as a function of 
the distance from the droplet center, hence c == 0. The 
nematic director q(°) is obtained from Eq.^ with R* — R 
The order parameter, defined in Eq.M, then simplifies to 



5(°)(r) 



(nj • q' 



(0) 



1. 



(13) 



^The constants in Eq.hd depend on the dimensionality of 
the system and are different from 3d, where, e.g. Q'"' — 

(3/2) J2i Hi (g) Hi - 1/2 holds. ~ 



The vector d is closely related to q'2) by a rotation oper- 
ation, where the rotation angle is the angle between q'2) 
and the x-axis. The direction d is where the field lines 
are radial; see the arrow in Fig.H. 



B. Defect distributions 



For a given configuration of particles the planar nematic 
droplet has a preferred direction given by the global ne- 
matic director q*^°) . Each of the two topological defects 
has a position c^ and an orientation di,z = 1,2. These 
quantities can be set in relation to each other to extract 
information about the average defect behavior and its 
fluctuations. In particular, we investigated the following 
probability distributions depending on a single distance 
or angle. 



Concerning single defect properties, we investigate the 
separation distance from the droplet center, 



P(r) = (2vrr)-ii^ (5(|c, 



r)), 



(18) 



i=l,2 



and the orientation relative to the nematic director, 

1=1,2 

Between both defects there is a distance distribution, 

P(ci2) = (27rci2)-i (5(|ci - C2I - C12)) , (20) 

and an angular distribution between defect orientations, 
P{0i2) = {S (arccos (di • d2) - ^12)) , (21) 



(-) (-) 
which can equivalently be defined with q^^^ ,q2 by us- 

(-) (-) 
ing the identity arccos(di ■ d2) = 2arccos(q]^^ • q2 )■ 



IV. COMPUTER SIMULATION 



A. Monte Carlo 



All our simulations were performed with the canonical 
Monte-Carlo technique keeping particle-number N, vol- 
ume V and temperature T constant, for details we re- 
fer to Ref. p^. To simulate spherocylinders with only 
hard interactions, each Monte-Carlo trial is exclusively 
accepted when there is no overlap of any particles. One 
trial always consists of a small variation of position and 
orientation of one HSC. 

For the planar case the translation for the particle i is 
constructed by adding a small random displacement Ar^ 
to the vector r^ , similarly the rotation consists of adding 
a small random vector An^ to the direction n^ with An^ ■ 

To achieve an isotropic trial on the surface of the sphere, 
the rotation matrix M is applied simultaneously to the 
vectors r^ and n^. It is defined as 

1 — c + a^c 7s -|- a/3c —(3s + ajc 

M:=( -7s + /3ac 1 - c + /3c a + f3jc | (22) 

ps + jac —as + ^(3c 1 — c + j'^c 



vector specifying the rotation axis, A9 is a small random 
angle. With this method a simultaneous translation and 
rotation is warranted by keeping the vectors r^ and n^ 
normalized and perpendicularly oriented. 

The maximal variation in all cases is adjusted such that 
the probability of accepting a move is about fifty percent. 
The overlap criteria were checked by comparing the sec- 
ond virial coefficient of two- and three-dimensional HSC 
with simulation results, where the excluded volume of 
two HSC were calculated. Each of the runs (I)- (VII) was 
performed with 5 • 10^ trials per particle. One tenth of 
each run was discarded for equilibration. Especially the 
strongly fluctuating distance distribution between both 
defects, P(ci2), needs good statistics. All quantities were 
averaged over 25 partial runs, from which also error bars 
were calculated. 

An overview of the simulated systems is given in Tab.|. 
The systems (I)- (VII) are planar. System (I) is the refer- 
ence. To study finite-size effects, system (II) has half as 
many particles, and system (III) has twice as many par- 
ticles as (I). To investigate the dependence on the ther- 
modynamic parameters, system (IV) has a lower packing 
fraction rj, and system (V) has a higher one compared to 
system (I). The other thermodynamic parameter is the 
anisotropy, which is smaller for system (VI) and higher 
for system (VII) compared to the system (I). To keep the 
nematic phase stable for the short rods of system (VI), 
the packing fraction 77 had to be increased. The pack- 
ing fraction of the dense system (V) is 77 = 0.4143. The 
spherical system has the same number of particles N, 
packing fraction rj and anisotropy p as the reference (I) . 
The radius of the sphere is half the radius of the planar 
droplet. The aligned rod model has the same parameters 
as the reference system (I). 



B. Technical issues 



We discuss briefly a projection method for the spherical 
problem and a search algorithm to find defect positions. 

In order to perform calculations for the spherical system 
all interesting vectors in three dimensions are projected 
to a two-dimensional plane. Imagine a given vector c 
from the middle of the sphere pointing to an arbitrary 
point of the surface. We convert a position Tj and orien- 
tation Ui to the vectors r? and n^ in a plane perpendic- 
ular to c through 



r^ = Ti - (C • TijC, 

nf = n, - (c • ni)c. 



(23) 
(24) 



with s — sin A6 and c — 1 — cos A9. a, (3, 7 are for every 
trial randomly chosen cartesian coordinates of the unit 



After obtaining a set {r^,n^} of three dimensional vec- 
tors on this way, we transform them into a set of two 



dimensional vectors by typical algebraic methods. As 
reference the projection of the x unit vector of the fixed 
three dimensional coordinate system is always the x- 
orientation of the "new" coordinate-system in two di- 
mensions. The results show that curvature effects are 
small. 

To investigate the radial structure and interactions of the 
disclinations it is necessary to localize the centers of the 
two point defects. As described in the last section, the 
A*- 2 '-parameter measures the degree of order of a half- 
integer defect in a chosen area, so the task is to find the 
two maxima of A'a-* in the droplet. In the planar case, we 
do this search with the following algorithm: A circular 
test-probe samples the droplet on a grid with a grid spac- 
ing of 5(7. At this points all the particles in the circle are 
taken to calculate A'- 2) in the described way. After sam- 
pling the grid both maxima are stored and for every max- 
imum a refining Monte-Carlo-search is performed. The 
surrounding of size of the grid spacing is randomly sam- 
pled and the probe is only moved when A' 2) increases. 
The search is stopped when the probe does not move for 
200 trials. In the spherical case the method is the same, 
but the grid is projected onto the sphere surface and the 
calculations of A^a' were performed with projected two- 
dimensional vectors as described before. 

It is important to chose an adequate radius R* for the 
probe. If R* is too large, the probe overlaps both de- 
fects. As they have opposite orientations on the aver- 
age, the located point of the maximum deviates from the 
point we are interested in. If the R* is too small, an ill- 
defined position results, as fluctuations become more im- 
portant. The simulation results show that a good choice 
is i?* = 2L. Although this definition contains some free- 
dom, we find the defect position to be a robust quantity. 
A detailed discussion is given in the following section. 



V. RESULTS 



A. Order within the droplet 



Let us discuss the order parameters S*'*' as a function 
of the radial distance from the center of the droplet; see 
Fig.^. S''^"-' is the usual bulk nematic order parameter, 
but radially resolved. It reaches values of 0.6-0.75 in the 
middle of the droplet, r < 2L, indicating a nematic por- 
tion that breaks the global rotational symmetry of the 
system. For r > 3L, 5'^°) decays to values slightly larger 
than the isotropic value of 0. The decrease, however, is 
not due to a microscopically isotropic fiuid state, as can 
be seen from the behavior of S^^\ This quantity indi- 
cates globally star-like alignment of particles for r > 3L. 
It vanishes in the nematic "street" in the center of the 
droplet. The distance where 5"*^°' and 5^^' intersect is 



an estimate for the defect positions. In Fig.^ the fi- 
nite size behavior of 5^*' is plotted for particle numbers 
N = 1004,2008,4016 corresponding to systems (II), (I), 
(III). There is a systematic shift of the intersection point 
of S*^"' and S'^^^ to larger values as the system grows, the 
numerical values are r/L — 2.54, 2.91, 3.87. However, if r 
is scaled by the droplet radius R, a slight shift to smaller 
values is observed as the system size grows. Keeping the 
medium-sized system (I) as a reference, we have inves- 
tigated the impact of changing the thermodynamic vari- 
ables. For different packing fractions, 77=0.2894 (IV), 
0.3321 (I), 0.4143 (V), we found that the intersection 
distances are r/L=3.90, 2.91, 1.43. In the bulk, upon in- 
creasing the density the nematic order grows. Here, this 
happens for the star-order S^^' . But this increase hap- 
pens on the cost of the nematic street (see S*^"-*) at small 
r- values. Increasing rj leads to a compression of the inho- 
mogeneous, interesting region in the center of the droplet. 
A similar effect can be observed upon changing the other 
thermodynamic variable, namely the anisotropy p. The 
nematic street is compressed for longer rods, p = 31 (VII), 
r/L=1.33. Shorter rods, p = 16, need a higher density 
to form a nematic phase, so the values for systems (I), 
r/L=3.16, and (VI), r/L=2.91, are similar, as both ef- 
fects cancel out. 

The behavior of S"*^^' is similar to the findings for a three- 
dimensional droplet, where a quadratic behavior near 
r = was predicted within Landau theory pl| ]. A sim- 
ulation study using the Lebwohl-Lasher model |3^ con- 
firmed this finding and revealed that a ring-like structure 
that breaks the spherical symmetry is present. A com- 
parison to the results for a 3d capillary by Andrienko and 
Allen [^ seems qualitatively possible as they find align- 
ment of particles predominantly normal to the cylinder 
axis. Their findings are consistent with the behavior of 
S^^' . Although our system is simpler as it only has two 
spatial dimensions, we could also establish the existence 
of a director field that breaks the spherical symmetry by 
considering the order parameter S^^' . 

Having demonstrated that the system exhibits a broken 
rotational symmetry, we have to assure that no freezing 
into a smectic or even crystalline state occurs. There- 
fore we plot radial density profiles p{r), where r is the 
distance from the droplet center, in Fig.O. The density 
shows pronounced oscillations for large r near the bound- 
ary of the system. They become damped upon increasing 
the separation distance from the droplet boundary and 
practically vanish after two rod lengths for intermediate 
density and four rod lengths for high density. Approach- 
ing the droplet center, r = 0, the density reaches a con- 
stant value for the weakly nematic systems (I) , (IV) , and 
(V). For the strongly nematic systems, (V) with high 
density and (VII) with large anisotropy, a density decay 
at the center of the droplet occurs. This effect is not di- 
rectly caused by the boundary as the density oscillations 
due to packing effects are damped. It is merely due to the 



topological defects present in the system. Quantitatively, 
the relative decrease is [p(3L) — p(0)] /p{3L) = 0.11 (V), 
0.09 (VII). The finite-size corrections for systems (II) and 
(III) are negligible. 

From both, the scissor-like behavior of the nematic or- 
der (Fig.^) and from the homogeneity of the density 
profile away from the system wall (Fig.g), we conclude 
that the system is in a thermodynamically stable nematic 
phase, and seems to contain two topological defects with 
charge 1/2. 

In a 2d bulk phase, two half-integer (1/2) defects are 
more stable than a single integer (1) defect, as the free 
energy is proportional to the square of the charge. How- 
ever, in the finite system of the computer simulation that 
is also affected by influence from the boundaries, it could 
also be possible that the defect pair merge into a single 
one 0,0. 

Next we investigate the defect positions and their orien- 
tations. To illustrate both, a snapshot of a configuration 
of the planar system is shown in Fig.0 (I). One can see 
the coupling of the nematic order from the first layer 
of particles near the wall to the inside of the droplet. 
The particles near the center of the droplet are aligned 
along a nematic director (indicated by the bar outside 
the droplet). The two emerging defects are depicted by 
symbols. See Fig.|| for a snapshot of the spherical sys- 
tem. There the total topological charge is not induced 
by a system boundary but by the topology of the sphere 
itself. 



B. Defect core 



The positions of the defects are defined by maxima of 
the A'-s) order parameter, see Section [II for its defini- 
tion. In Fig.^ A*-2 ^ is plotted as a function of the spatial 
coordinates r^ and Vy for one given configuration. There 
are two pronounced maxima, indicated by bright areas, 
which are identified as the positions of the defect cores 
Ci and C2. There are several more local maxima appear- 
ing as gray islands. These are identified as statistical 
fluctuations already present in the bulk nematic phase. 



A drift of the positions of a defect core was also reported 
in |32|. Here we follow this motion, to investigate the 
surrounding of the defects. The order parameter S'^^) 
is radially resolved around the defect position in FigjlC|. 
It has a pronounced maximum around r = 1.2L. For 
smaller distances it decreases rapidly due to disorder in 
the core region. For larger distances the influence from 
the second defect partner decreases the half-integer or- 
der S^'^K Increasing the overall density, and increasing 
the anisotropy leads to a more pronounced hump. The 
finite-size corrections, (II), (HI), and the boundary effects 



(sphere) are negligible. However, the curves show two ar- 
tifacts: A rise near r — and a jump at the boundary of 
the search probe, r = 2L. In the inset the profile around 
a bulk defect is shown. It has a plateau value inside the 
probe, r < 2L, and vanishes outside. If we subtract this 
contribution from the pure data (I) , continuous behavior 
at r = 2L can be enforced. 

However, the model does not account for 3d effects like 
the "biaxial escape" , namely the sequence planar uniaxial 
- biaxial - uniaxial with increasing distance from the core 
center |3J] , as the particles are only 2d rotators. Schopohl 
and Sluckin pO| found an interface-like behavior between 
the inner and outer parts of a disclination line in 3d. In 
our system, we do not find a sign of an interface between 
the isotropic core and the surrounding nematic phase. 
This might be due to a small interface tension and a 
very weak bulk nematic-isotropic phase transition. 

By radially resolving the probability of finding a particle 
around a defect center, we end up with density profiles 
depicted in Fig.O. The defect is surrounded by density 
oscillations with a wavelength of the particle length. The 
finite-size dependence is small. To estimate the influence 
from the system wall, one may compare with the spher- 
ical system. It shows slightly weaker oscillations. This 
might be due to curvature effects, as the effective packing 
fraction is slightly smaller as the linear particles may es- 
cape the spherical system. The toy model of aligned rods 
also exhibits a non-trivial density profile, showing a de- 
crease towards small distance and oscillations compared 
to rotating rods. In all cases the first peak has a sepa- 
ration distance of half a particle length from the defect 
center. The second peak appears at r = 3/2L. Again the 
search probe induces an artificial structure near r = 2L. 
From this analysis, we can conclude that the oscillations 
are due to packing effects. The density oscillations be- 
come more pronounced at higher density, and for larger 
anisotropy, see Fig.O. 



C. Defect position 



In the planar system, each defect is characterized by its 
radial distance r from the center, and the angle 9 be- 
tween its orientation and the global nematic director q(°) . 
We discuss the probability distributions of these quanti- 
ties. In Fig.^ the distribution for finding the defect at 
a distance r from the center is shown. Generally, the 
distributions are very broad. This indicates large mobil- 
ity of the defects. Changing the thermodynamical vari- 
ables has a large effect. For the stronger nematic sys- 
tems (V) and (VII), the distribution becomes sharper 
with a pronounced maximum at r = 1.5L. Decreas- 
ing the anisotropy weakens the nematic phase, so system 
(IV) has a very broad distribution. The inset shows that 



the distribution becomes broader upon increasing system 
size. 



E. Outlook 



D. Interactions between two defects 



A complete probability distribution of both positions of 
the defect cores can be regarded as arising from an effec- 
tive interaction potential T4ff(ci, C2) between the defects. 
The latter play the role of quasi-particles. The effective 
interaction arises from averaging over the particle posi- 
tions while keeping the defect positions constant. The 
effective interaction and the probability distribution are 
related via P(ci,C2) oc exp(— /3V^ff(ci, C2)). 

Instead of the full probability distribution, we show its 
dependence on the separation distance between both de- 
fects and on their relative orientation. In FigjlJ the prob- 
ability distribution of finding two defects at a distance 
C12 is shown. It has small values for small as well as 
large C12. Hence at small distances the defects repel each 
other. At large distances their effective interaction is at- 
tractive. Increasing the nematic order by increasing the 
density (V) or rod length (VII) causes the average defect 
separation distance to shrink. The rise near r / L = 1 is 
an artifact: These are events where the search algorithm 
does not find two different defects, but merely finds the 
same defect two times. To avoid the problem a cutoff at 
r = L was introduced. The finite size behavior is strong; 
see the inset. The large system (III) allows the defects 
to move further away from each other, whereas in the 
smaller system (II) they are forced to be closer together. 
However, from the simulation data, it is hard to obtain 
the behavior in the limit R/ L —>■ 00. 

This is somewhat in contrast to the phase diagram of a 
3d capillary pM containing isotropic, planar-radial and 
planar-polar structures, if one is willing to identify the 
dependence on temperature with our athermal system. 
There it was found that the transition from the planar- 
polar to the planar-radial structure happens upon in- 
creasing the temperature (and hence decreasing the ne- 
matic order). 

The difference angle 6*12 between both defect orientations 
in the planar system, see Fig.|l^, is most likely tt, hence 
the defects point on average away from each other. How- 
ever, the orientations are not very rigid. For the least 
ordered system (IV) there is still a finite probability of 
finding the defects with a relative orientation of 90 de- 
grees! Even for the strongly nematic systems (V) and 
(VII) the angular fluctuations are quite large. The inset 
in Fig.|l^ shows the distribution of the angle 9 between 
the defect orientation and the global nematic director. A 
clear maximum near tt/2 occurs. Again, the distributions 
become sharper as density or anisotropy increase. 



Finally, it is worth mentioning that the spherical sys- 
tem still contains surprises. See Fig. 16 for an unexpected 
configuration, namely an assembly of three positive 1/2- 
defects sitting at the corners of a triangle and a negative 
-1/2-defect in its center. This is remarkable, because the 
negative defect could annihilate with one of the outer 
positive defects. 

In all cases, integer defects seem to dissociate into half- 
integer defects. The complete equilibrium defect distri- 
bution of hard sphero cylinders lying tangentially on a 
sphere remains an open question. 



VI. CONCLUSIONS 



In conclusion, we have investigated the microscopic struc- 
ture of topological defects of nematics in a spherical 
droplet with appropriate homeotropic boundary and for 
particles lying on the surface of a sphere. We have used 
hard spherocylinders as a model system for a lyotropic 
nematic liquid crystal. This system allows us to study 
the statistical behavior of the microscopic rotational and 
positional degrees of freedom. For this system we find 
half-integer topological point defects in two dimensions 
to be stable. The defect core has a radius of the order 
of one particle length. As an important observation, the 
defect generates a free-standing density oscillation. It 
possesses a wavelength of one particle length. Consid- 
ering the defects as fluctuating quasi-particles we have 
presented results for their effective interaction. 

The microscopic structure revealed by radially resolving 
density and order parameter profiles around the defect 
position is identical for the planar and the spherical sys- 
tem. 

An experimental investigation using anisotropic colloidal 
particles fBjjG^l like tobacco mosaic viruses or carbon 
nanotubes is highly desirable to test our theoretical pre- 
dictions. Then larger accessible system sizes can be ex- 
ploited. Also of interest is the long-time dynamical be- 
havior of the motion of topological defects. The advan- 
tage of colloidal systems over molecular liquid crystals is 
the larger length scale that enables real-space techniques 
like digital video-microscopy to be used. 

From a more theoretical point of view it would be inter- 
esting to describe the microstructure of topological de- 
fects within the framework of density functional theory. 
Using phenomenological Ginzburg-Landau models, one 
could take the elastic constants of the HSC model as an 
input, and could calculate the defect positions and check 
against our simulations. 



Finally we note that we currently investigate the three- 
dimensional droplets that are filled with spherocylinders. 
In this case more involved questions appear, as both, 
point and line defects, may appear. 
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FIG. 1. Two hard spherocylinders with position coordi- 
nates Ti and Tj , and orientations n^ and rij . The width of the 
particles is a, the total rod length is denoted by L. 
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TABLE I. Overview of the simulated parameter range: 
number of particles A'', anisotropy p, packing fraction r), scaled 
droplet diameter 2R/ L. Systems (I)-(VII) are planar, the 
system named "sphere" corresponds to spherical geometry. 



FIG. 2. Homeotropic boundary conditions for the planar 
droplet. The particle centers (points) are not allowed to cross 
a circle with diameter R — L/2 (dashed line) . Then the shape 
of each particle lies inside a circle with radius R. 




FIG. 3. Spherical system. Each particle with position r^ 
and orientation n^ is forced to lie tangentially on the surface 
of a sphere. 
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FIG. 4. Model of aligned rods. Each particle (discorect- 
angles) has an orientation according to a prescribed director 
field (lines). The position of the arising 1/2-defect is indicated 
by a filled circle, the orientation by an arrow. 




FIG. 6. Radially resolved density profiles p{r) as a func- 
tion of the distance from the droplet center r scaled by the 
particle length L. System (I) is reference, compared to lower 
(IV) and higher (V) packing fraction, and lower (VI) and 
higher (VII) anisotropies. The inset shows the behavior near 
the origin where a density decrease for (V) and (VII) appears 
for r < 21/. 
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FIG. 7. Snapshot of a typical particle configuration for 
the planar system (I). The particles are rendered dark. The 
two black symbols inside the droplet indicate positions and 
orientations of defects. The black bar outside the droplet 
indicates the global nematic director q'°'. 



FIG. 5. Nematic order parameters S'*' as a function of 
the radial distance r from the droplet center, scaled by the 
rod length L. Star order 5' ' and bulk order S^ ' is shown. 
System (I) is reference, (II) has halved and (III) has dou- 
bled particle number. See Tab.| for a compilation of system 
parameters. Error bars are only given for (I). 




FIG. 8. Snapshot of a typical particle configuration for 
the spherical system. The particles are rendered dark. There 
is one 1/2-defect on the left side and one on the right side. 
They point away from each other. 
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FIG. 9. Order parameter A'2) as a function of spatial 
coordinates rx,ry. Bright areas correspond to large values, 
dark areas correspond to small values of A' 2 ) . The two bright 
spots near the center are identified as topological defects, the 
gray islands as bulk defects. 



FIG. 11. Density profile as a function of the distance 
from the defect center. System (I) is reference, (II) has fewer 
particles, (III) has more. The spherical and aligned models 
are shown. 
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FIG. 10. Order parameter profiles S^'^^ around the defect 
center as a function of the scaled distance r/L from the de- 
fect center. The reference system (I) is to be compared with 
lower (IV) and higher (V) packing fraction, and lower (VI) 
and higher (VII) anisotropics. The inset shows S^^' for bulk 
defects and for the difference between (I) and the bulk. 
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FIG. 12. Same as Figjll], but for lower (IV) and higher 
(V) packing fraction and shorter (VI) and longer particles 
(VII), compared to system (I). 
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FIG. 13. Probability distribution P{r) for the distance of 
a defect from the center of the droplet r/L, for lower (IV) and 
higher (V) packing fraction and shorter (VI) and longer par- 
ticles (VII), compared to system (I). The inset shows the the 
finite size behavior for halved (II) and doubled (III) particle 
number. 
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FIG. 16. Triangular configuration of three positive de- 
fects around a spontaneously formed negatively charged de- 
fect (central dot). 



FIG. 14. Probability distribution P(ci2) for the separa- 
tion distance between both defect positions scaled by the par- 
ticle length for lower (IV) and higher (V) packing fraction and 
shorter (VI) and longer spherocylinders (VII) , as compared to 
system (I). The inset shows the finite size behavior for halved 
(II) and doubled (III) particle number compared to (I). 




FIG. 15. Probability distribution P{di2) for the difference 
angle between both defect orientations. The reference system 
(I) is to be compared with lower (IV) and higher (V) pack- 
ing fraction, and lower (VI) and higher (VII) anisotropies. 
The inset shows the distribution P{0) of the difference an- 
gle between the direction of one of the defects and the global 
nematic director for the same parameters. 
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